doubleprecision Function omegalaw(r) result(chi)
!     r - cylinder radius, r < 1 inside the star
!     w=1/(1+mult*r^n),
!     chi=rot. potential
 Implicit none
 real(8) r,x1
 integer, parameter :: n=4 ! the power. n must be 1, 2, or 4
 real(8), parameter :: ratio=10.d0 ! The ratio=w(0)/w(1)
 real(8), parameter :: mult=ratio-1
!--------------------------------------------------------------------

 select case(n)
  case(1) !--------------/1/------------------------
   x1=mult*dmin1(r,1.d0)
   chi=(dlog(1+x1)-x1/(1+x1))/mult**2

  case(2) !-------------/2/------------------------
   x1=mult*dmin1(r**2,1.d0)
   chi=0.5d0*x1/(1+x1)/mult

  case(4) !-------------/4/------------------------
   x1=dsqrt(mult)
   chi=0.25d0/x1
   x1=x1*dmin1(r**2,1.d0)
   chi=chi*(datan(x1)+x1/(1+x1**2))

  case default !-------/exit with wrong value of n/--------------
   print*,'wrong values if omgalaw, n=',n
   read*
   stop
 end select

 return

end function

!=========================================================================
!*************************************************************************

doubleprecision Function savedlaw(r) result(chi)
!     r - cylinder radius, r < 1 inside the star
!
 Implicit none
 integer,parameter :: ndot=150
 real(8) r,r2,x1,x2,x3
 real(8),save :: rs(ndot),rs2(ndot),w2(ndot),intg(ndot)
 integer, save :: first
 integer k,i

 data first/0/
! data rs/0.00666667d0,0.0133333d0,0.02d0,0.0266667d0,0.0333333d0,0.04d0,&
!         0.0466667d0,0.0533333d0,0.06d0,0.0666667d0,0.0733333d0,0.08d0,&
!         0.0866667d0,0.0933333d0,0.1d0,0.10667d0,0.11333d0,0.12d0,0.12667d0,&
!         0.13333d0,0.14d0,0.14667d0,0.15333d0,0.16d0,0.16667d0,0.17333d0,&
!         0.18d0,0.18667d0,0.19333d0,0.2d0,0.20667d0,0.21333d0,0.22d0,&
!         0.22667d0,0.23333d0,0.24d0,0.24667d0,0.25333d0,0.26d0,0.26667d0,&
!         0.27333d0,0.28d0,0.28667d0,0.29333d0,0.3d0,0.30667d0,0.31333d0,&
!         0.32d0,0.32667d0,0.33333d0,0.34d0,0.34667d0,0.35333d0,0.36d0,0.36667d0,&
!         0.37333d0,0.38d0,0.38667d0,0.39333d0,0.4d0,0.40667d0,0.41333d0,0.42d0,&
!		 0.42667d0,0.43333d0,0.44d0,0.44667d0,0.45333d0,0.46d0,0.46667d0,&
!		 0.47333d0,0.48d0,0.48667d0,0.49333d0,0.5d0,0.50667d0,0.51333d0,&
!		 0.52d0,0.52667d0,0.53333d0,0.54d0,0.54667d0,0.55333d0,0.56d0,0.56667d0,&
!		 0.57333d0,0.58d0,0.58667d0,0.59333d0,0.6d0,0.60667d0,0.61333d0,0.62d0,&
!		 0.62667d0,0.63333d0,0.64d0,0.64667d0,0.65333d0,0.66d0,0.66667d0,&
!		 0.67333d0,0.68d0,0.68667d0,0.69333d0,0.7d0,0.70667d0,0.71333d0,&
!		 0.72d0,0.72667d0,0.73333d0,0.74d0,0.74667d0,0.75333d0,0.76d0,&
!		 0.76667d0,0.77333d0,0.78d0,0.78667d0,0.79333d0,0.8d0,0.80667d0,&
!		 0.81333d0,0.82d0,0.82667d0,0.83333d0,0.84d0,0.84667d0,0.85333d0,&
!		 0.86d0,0.86667d0,0.87333d0,0.88d0,0.88667d0,0.89333d0,0.9d0,0.90667d0,&
!		 0.91333d0,0.92d0,0.92667d0,0.93333d0,0.94d0,0.94667d0,0.95333d0,&
!		 0.96d0,0.96667d0,0.97333d0,0.98d0,0.98667d0,0.99333d0,1.d0/
! data w2/0.56567d0,0.64117d0,0.72074d0,0.75074d0,0.76422d0,0.77084d0,0.77413d0,&
!         0.77558d0,0.77594d0,0.77559d0,0.77472d0,0.77349d0,0.77196d0,0.77018d0,&
!		 0.76819d0,0.76601d0,0.76365d0,0.76113d0,0.75845d0,0.75563d0,0.75266d0,&
!		 0.74956d0,0.74632d0,0.74294d0,0.73944d0,0.73581d0,0.73206d0,0.72819d0,&
!		 0.72419d0,0.72009d0,0.71586d0,0.71153d0,0.70709d0,0.70254d0,0.69789d0,&
!		 0.69314d0,0.68829d0,0.68334d0,0.6783d0,0.67317d0,0.66796d0,0.66266d0,&
!		 0.65727d0,0.65181d0,0.64627d0,0.64066d0,0.63498d0,0.62923d0,0.62342d0,&
!		 0.61754d0,0.6116d0,0.60561d0,0.59957d0,0.59347d0,0.58733d0,0.58114d0,&
!		 0.57491d0,0.56865d0,0.56234d0,0.556d0,0.54964d0,0.54324d0,0.53682d0,&
!		 0.53038d0,0.52391d0,0.51743d0,0.51094d0,0.50444d0,0.49792d0,0.4914d0,&
!		 0.48487d0,0.47835d0,0.47182d0,0.4653d0,0.45878d0,0.45227d0,0.44577d0,&
!		 0.43928d0,0.4328d0,0.42635d0,0.41991d0,0.41349d0,0.40709d0,0.40072d0,&
!		 0.39438d0,0.38806d0,0.38177d0,0.37552d0,0.3693d0,0.36311d0,0.35696d0,&
!		 0.35085d0,0.34478d0,0.33874d0,0.33276d0,0.32681d0,0.32091d0,0.31506d0,&
!		 0.30925d0,0.3035d0,0.29779d0,0.29214d0,0.28653d0,0.28098d0,0.27549d0,&
!		 0.27005d0,0.26466d0,0.25933d0,0.25406d0,0.24885d0,0.2437d0,0.2386d0,&
!		 0.23357d0,0.2286d0,0.22369d0,0.21884d0,0.21405d0,0.20933d0,0.20467d0,&
!		 0.20007d0,0.19554d0,0.19108d0,0.18668d0,0.18234d0,0.17807d0,0.17387d0,&
!		 0.16974d0,0.16567d0,0.16167d0,0.15773d0,0.15387d0,0.15007d0,0.14634d0,&
!		 0.14268d0,0.13909d0,0.13557d0,0.13211d0,0.12873d0,0.12542d0,0.12218d0,&
!		 0.11901d0,0.11591d0,0.11289d0,0.10994d0,0.10706d0,0.10426d0,0.10154d0,&
!		 0.0988922d0,0.0963299d0,0.0944775d0/

!--------------------------------------------------------------------

 if(first.eq.0)then ! first call
  first=1

!------/reading/-----------

 open(303,file='Omega.dat',err=11)
 do k=1,ndot
  read(303,*,err=11)i,rs(k),x1,x2,x3,w2(k)
 end do
 close(303)

!--------------------------

  rs2=rs**2
  intg(1)=w2(1)*rs2(1)
  do k=2,ndot
   intg(k)=intg(k-1)+0.5d0*(rs2(k)-rs2(k-1))*(w2(k)+w2(k-1))
  end do
  intg=0.5d0*intg
 end if

 if(r.le.rs(1))then
  chi=0.5d0*w2(1)*r**2
  return
 elseif(r.ge.1.d0)then
  chi=intg(ndot)
  return
 end if

 k=minloc(r-rs,mask=r-rs.ge.0.d0,dim=1)

 r2=r**2
 chi=intg(k)+0.25d0*(r2-rs2(k))/(rs2(k+1)-rs2(k))*&
     (w2(k)*(2*rs2(k+1)-rs2(k)-r2)+w2(k+1)*(r2-rs2(k)))

 return

11 print*,'Error with data-file Omega.dat'
read*
stop

end function

